9. Γεωεπεξεργασία διανυσματικών δεδομένων

Ο χώρος και οι γεωμετρικές δομές του μπορούν να αναπαρασταθούν μέσω διανυσματικών δεδομένων (Vector). Οι τρείς βασικές δομές δεδομένων είναι:

  • Τα σημεία

  • Οι γραμμές

  • Τα πολύγωνα

alt text

Επιπλέον αυτά τα γεωμετρικά δεδομένα συνοδεύονται και απο περιγραφικά δεδομένα που αφορούν τις ιδιότητες ή τα χαρακτηριστικά αυτών των δεδομένων.

alt text

Στα Σ.Γ.Π. ο πιο συνηθισμένος τύπος αρχείων αποθήκευσης αυτών των δεδομένων είναι το shapefile. Πλέον έχουν αναπτυχθεί και άλλοι τύποι όπως geojson, geopackage και χωρικές βάσεις (geodatabases) όπως η Postgresql/Postgis.

Ο προσδιορισμός της γεωγραφικής θέσης στην υδρόγειο γίνεται μέσω ενός ζεύγους γεωγραφικών συντεταγμένων. Κάθε σημείο στον χώρο προσδιορίζεται γεωγραφικά από το γεωγραφικό μήκος (λ) και το γεωγραφικό πλάτος (φ).

alt text

Για να αποδοθεί η τρισδιάστατη υδρόγειος σφαίρα σε ένα δυσδιάστατο σύστημα αναφορά χρησιμοποιείται ένα προβολικό σύστημα.

alt text

Κάθε προβολικό σύστημα της σφαίρας στο επίπεδο εισάγει μια σειρά παραμορφώσεων που αφορά το σχήμα των γεωμετρικών δομών, την κλίμακα, την έκταση και τις αποστάσεις. Ανάλογα το προβολικό σύστημα κάποιες από τις παραπάνω παραμορφώσεις εμφανίζονται σε μεγάλο βαθμό και άλλες όχι. Οπότε ανάλογα το είδος της έρευνας ο ερευνητής οφείλει να γνωρίζειπαραμορφώσεις εισάγει η κάθε προβολή και ανάλογα να επιλέγει την προβολή με τα λιγότερα σφάλματα.

Python βιβλιοθήκες για διανυσματικά δεδομένα

Για την ανάγνωση, εγγραφή, επεξεργασία διανυσματικών δεδομένων στην Python έχουν καθιερωθεί μια σειρά βιβλιοθηκών. Η αρχαιότερη και βασική βιβλιοθήκη είναι η GDAL/OGR. Επειδή η βιβλιοθήκη δεν είναι ιδιαίτερα συμβατή σε σχέση με τον τρόπο συγγραφής της Python και επιρρεπής σε σφάλματα έχουν αναπτυχθεί πιο σύγχρονες βιβλιοθήκες όπως η βιβλιοθήκη [Fiona] που είναι ιδιαίτερα χρήσιμη για την ανάγνωση/εγγραφή διανυσματικών δεδομένων Fiona και η βιβλιοθήκη Shapely η οποία χρησιμοποιείται για την επεξεργασία και ανάλυση.

Η βιβλιοθήκη Shapely

shapely from WKT

Μπορούμε να δημιουργήσουμε αντικείμενα shapely που αναπαριστούν σημεία ή γραμμές ή πολύγωνα μέσω WKT. Η γλώσσα (WKT)[https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry] είναι μια ειδική διάλεκτος για την περιγραφή διανυσματικών αντικειμένων. Εισάγουμε τις απαραίτητες υπο-βιβλιοθήκες (geometry και wkt) από την βιβλιοθήκη shapely.

import shapely.geometry
import shapely.wkt

Πολύγωνα

Καλούμε την μέθοδο shapely.wkt.loads() για να δημιουργήσουμε shapely objects από wkt

pol1 = shapely.wkt.loads("POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0,0 0))")
pol1
../_images/Lesson_9_19_0.svg

Η εκτύπωση του shapely αντικειμένου επιστρέφει μία περιγραφή σε μορφή WKT.

print(pol1)
POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0))

Εκτύπωση του τύπου του αντικειμένου pol1

type(pol1)
shapely.geometry.polygon.Polygon
pol2 = shapely.wkt.loads("POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))")
print(pol2)
POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))
pol2
../_images/Lesson_9_25_0.svg

Δημιουργία MultiPolygon shapely object από αντίστοιχη συλλογή πολυγώνων MULTIPOLYGON

pol3 = shapely.wkt.loads("""
MULTIPOLYGON
(((40 40, 20 45, 45 30, 40 40)),
((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))
""")
type(pol3)
shapely.geometry.multipolygon.MultiPolygon
pol3
../_images/Lesson_9_29_0.svg
type(pol3)
shapely.geometry.multipolygon.MultiPolygon

Σημεία

pnt1 = shapely.wkt.loads("""
POINT (30 10)
""")
pnt1
../_images/Lesson_9_33_0.svg

Συλλογή σημείων

pnt2 = shapely.wkt.loads("""
MULTIPOINT(0 0,1 1)
""")
pnt2
../_images/Lesson_9_36_0.svg

Γραμμές

line1 =shapely.wkt.loads("""
LINESTRING(1.5 2.45,3.21 4)
""")

line1
../_images/Lesson_9_38_0.svg

Συλλογή γραμμών

line2 =shapely.wkt.loads("""
MULTILINESTRING((0 0,-1 -2,-3 -4, -3 -8),(2 3,3 4,6 7))
""")

line2
../_images/Lesson_9_40_0.svg

Shapely μέσω συναρτήσεων

from shapely.geometry import Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon

Κάθε συνάρτηση λαμβάνει διαφορετικά ορίσματα ανάλογα με την συνάρτηση.

Στην συνέχεια ακολουθεί η δημιουργία ενός αντικειμένου shapely γεωμετρία σημείου.

pnt1 = shapely.geometry.Point((2, 0.5))
pnt1
../_images/Lesson_9_44_0.svg

Για την δημιουργία συλλογής σημείων multipoint δίνουμε σαν όρισμα στην συνάρτηση shapely.geometry.MultiPoint μία λίστα από πλειάδες (tuples). Η κάθε πλειάδα (x,y)αντιστοιχεί στις συντεταγμένες ενός σημείου.

coords = [(2, 0.5), (1, 1), (-1, 0), (1, 0)]
pnt2 = shapely.geometry.MultiPoint(coords)
pnt2
../_images/Lesson_9_46_0.svg

Για την δημιουργία γραμμών χρησιμοποιείται πάλι μία λίστα από tuples που περιγράφουν τις κορυφές της γραμμής. Στο παρακάτω παράδειγμα χρησιμοποιούμαι την προηγούμενη πλειάδα αλλά πλέον χρησιμοποιούμε την μέθοδο shapely.geometry.LineString για την δημιουργία γραμμής και όχι την μέθοδο hapely.geometry.MultiPoint που δημιουργεί συλλογές σημείων.

line1 = shapely.geometry.LineString(coords)
line1
../_images/Lesson_9_48_0.svg

Κατασκεύη γραμμής από Point Objects:

p1 = shapely.wkt.loads("POINT (0 0)")
p2 = shapely.wkt.loads("POINT (1 1)")
p3 = shapely.wkt.loads("POINT (2 -1)")
p4 = shapely.wkt.loads("POINT (2.5 2)")
p5 = shapely.wkt.loads("POINT (1 -1)")

shapely.geometry.LineString([p1, p2, p3, p4, p5])
../_images/Lesson_9_50_0.svg

Αντίστοιχα μπορούμε να δημιουργήσουμε συλλογές γραμμών.

Δημιουργούμε ξεχωριστά αντικείμενα γραμμών shapely και στην συνέχεια καλούμε την συνάρτηση shapely.geometry.MultiLineString όπου ορίζουμε σαν όρισμα μια λίστα με τις μεμονωμένες γραμμές.

l1 = shapely.geometry.LineString([(2, 0.5), (1, 1), (-1, 0), (1, -1)])
l2 = shapely.geometry.LineString([(-2, 1), (2, -0.5), (3, 1)])
line2 = shapely.geometry.MultiLineString([l1, l2])
line2
../_images/Lesson_9_52_0.svg

Αντίστοιχα δημιουργούμε πολυγωνικά αντικείμενα μέσω μια λίστας με πλειάδες σημείων που χρησιμοποιείται σαν όρισμα στην συνάρτηση shapely.geometry.Polygon

coords = [(0, 0), (0, -1), (7.5, -1), [7.5, 0], (0, 0)]
shapely.geometry.Polygon(coords)
../_images/Lesson_9_54_0.svg

Αν θέλουμε μπορούμε να περάσουμε μία δεύτερη λίστα η οποία περιλαμβανει επιμέρους λίστες με πλειάδες σημείων τα οποία περιγράφουν τρύπες μέσα στο πολύγωνο.

coords_exterior = [(0, 0), (0, -1), (7.5, -1), (7.5, 0), (0, 0)]
coords_interiors = [[(0.4, -0.3), (5, -0.3), (5, -0.7), (0.4, -0.7), (0.4, -0.7)]]
shapely.geometry.Polygon(coords_exterior, coords_interiors)
../_images/Lesson_9_56_0.svg

Με την συνάρτηση shapely.geometry.MultiPolygon δημιουργούμε αντίστοιχα συλλογές πολυγώνων από μεμονωμένα αντικείμενα shapely.geometry.polygon.Polygon

multipolygon1 = shapely.geometry.MultiPolygon([pol2, pol3])
multipolygon1
/home/leonidas/anaconda3/envs/book/lib/python3.10/site-packages/shapely/geometry/multipolygon.py:202: ShapelyDeprecationWarning: __getitem__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  shell = ob[0]
/home/leonidas/anaconda3/envs/book/lib/python3.10/site-packages/shapely/geometry/multipolygon.py:203: ShapelyDeprecationWarning: __getitem__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.
  holes = ob[1]
../_images/Lesson_9_58_1.svg
print(multipolygon1)
MULTIPOLYGON (((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1)), ((40 40, 20 45, 45 30, 40 40)))

Σύμφωνα με τις προδιαγραφές simple features το παραπάνω object δεν είναι έγκυρο γιατί ένα πολύγωνο τέμνει ένα άλλο σε άπειρο αριθμό σημείων. Μπορούμε να ελέγξουμε την εγκυρότητα ενός αντικειμένου με την κλήση τις ιδιότητας is_valid. Επίσης κατά την οπτικοποίηση στο προηγούμενο βήμα του αντικειμένου multipolygon1 αυτό εμφανίζεται με κόκκινο (αντί για πράσινο). Ένδειξη ότι δεν ειναι valid.

multipolygon1.is_valid
True

Ταυτόχρονα μπορούμε να φτιάξουμε σύνθετες συλλογές από επιμέρους αντικείμενα shapely.

geo_collection = shapely.geometry.GeometryCollection([multipolygon1,line2,pnt1])
geo_collection
../_images/Lesson_9_64_0.svg

Shapely μέσω shape

Μπορούμε να δημιουργήσουμε αντικείμενα shapely μέσω της συνάρτησης shapely.geometry.shape η οποία δέχεται σαν όρισμα ένα λεξικό μορφής GEOJSON το οποίο πρέπει να έχει τις εξής δύο ιδιότητες.

  • Την ιδιότητα "type" που περιγράφει τον τύπο γεωμετρίας

  • Την ιδιότητα "coordinates" που περιγράφει τις γεωμετρίες και τις συντεταγμένες τους σαν λίστες ή πλειάδες.

d = {"type": "Point", "coordinates": (0, 1)}
shapely.geometry.shape(d)
../_images/Lesson_9_67_0.svg

Δημιουργία αντικειμένου MultiPolygon μέσω GEOJSON λεξικού:

d = {
  "type": "MultiPolygon",
  "coordinates": [
    [
      [[40, 40], [20, 45], [45, 30], [40, 40]]
    ],
    [
      [[20, 35], [10, 30], [10, 10], [30, 5], [45, 20], 
      [20, 35]], 
      [[30, 20], [20, 15], [20, 25], [30, 20]]
    ]
  ]
}
pol3 = shapely.geometry.shape(d)
pol3
../_images/Lesson_9_69_0.svg

Γεωμετρικός τύπος δεδομένων

H ιδιότητα .geom_type ενός αντικειμένου shapely περιγράφει τον γεωμετρικό τύπο της:

pol1.geom_type
'Polygon'
line1.geom_type
'LineString'
geo_collection.geom_type
'GeometryCollection'

Συντεταγμένες

Για να ανακτήσουμε τις συντεταγμένες ενός shapely αντικειμένου καλούμε με διαφορετικό τρόπο τις ανάλογες συναρτήσεις ανάλογα την πολυπλοκότητα του κάθε τύπου.

Για ένα Point object καλούμε άμεσα την ιδιότητα coords η οποία επιστρέφει ένα shapely.coords.CoordinateSequence object.

pnt1.coords
<shapely.coords.CoordinateSequence at 0x7f970f4f2ad0>

Μπορούμε να λάβουμε ως λίστα τις πλειάδες συντεταγμένων που το απαρτίζουν

list(pnt1.coords)
[(2.0, 0.5)]

Αντίστοιχα για γραμμή

list(line1.coords)
[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, 0.0)]

Για να ελέγξουμε σε ένα αντικείμενο συλλογών γεωμετρίας (MultiPoint, MultiPolygon κτλ) πόσες επιμέρους γεωμετρίες περιλαμβάνει:

len(line2.geoms)
2
line2
../_images/Lesson_9_85_0.svg

Για να λάβουμε το πρώτο αντικείμενο γεωμετρίας:

line2.geoms[0]
../_images/Lesson_9_87_0.svg
type(line2.geoms[0])
shapely.geometry.linestring.LineString
line2.geoms[0].geom_type
'LineString'

Και λαμβάνουμε τις συντεταγμένες της πρώτης γεωμετρίας:

list(line2.geoms[0].coords)
[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, -1.0)]

ή της δεύτερης

list(line2.geoms[1].coords)
[(-2.0, 1.0), (2.0, -0.5), (3.0, 1.0)]

Για τα πολύγωνα ακολουθούμε διαφορετική προσέγγιση. Ένα πολύγωνο αποτελεί από το εξωτερικό περίγραμμα (exterior) ή και ένα ή περισσότερα περιγράμματα εσωτερικών τρυπών (interiors). Κατά συνέπεια έχουμε συντεταγμένες που περιγράφουν το κάθε περίγραμμα.

Το εξωτερικό περίγραμμα του pol1 αντικειμένου:

pol1.exterior
../_images/Lesson_9_96_0.svg

Και οι συντεταγμένες του:

pol1
../_images/Lesson_9_98_0.svg
list(pol1.exterior.coords)
[(0.0, 0.0), (0.0, -1.0), (7.5, -1.0), (7.5, 0.0), (0.0, 0.0)]

Το pol1 όμως δεν έχει εσωτερικές τρύπες γι αυτό και το παρακάτω επιστρέφει μηδέν.

len(pol1.interiors)
0

Έστω το παρακάτω MultiPolygon object που δημιουργήσαμε σε προηγούμενο στάδιο.

pol3
../_images/Lesson_9_103_0.svg

Περιλαμβάνει δύο ξεχωριστά γεωμετρικά αντικείμενα:

len(pol3.geoms)
2

Ας πάρουμε το πρώτο:

pol3.geoms[0]
../_images/Lesson_9_107_0.svg

Δεν περιλαμβάνει καμία τρύπα στο εσωτερικό του. Ας πάρουμε τις συντεταγμένες από το περίγραμμά του (exterior)

pol3.geoms[0].exterior.coords
<shapely.coords.CoordinateSequence at 0x7f9750d313c0>

Και ας τις επιστρέψουμε σαν λίστα πλειάδων από ζεύγη της μορφής (x,y)

list(pol3.geoms[0].exterior.coords)
[(40.0, 40.0), (20.0, 45.0), (45.0, 30.0), (40.0, 40.0)]

Ας δοκιμάσουμε το δεύτερο αντικείμενο γεωμετρίας. Ας το δούμε:

pol3.geoms[1]
../_images/Lesson_9_113_0.svg

Λαμβάνουμε τις συντεταγμένες του εξωτερικού περιγράμματος:

list(pol3.geoms[1].exterior.coords)
[(20.0, 35.0),
 (10.0, 30.0),
 (10.0, 10.0),
 (30.0, 5.0),
 (45.0, 20.0),
 (20.0, 35.0)]

Ας δούμε πόσες τρύπες έχει στο εσωτερικό του:

len(pol3.geoms[1].interiors)
1

Παίρνουμε το περίγραμμα της τρύπας

pol3.geoms[1].interiors[0]
../_images/Lesson_9_119_0.svg

Και τις συντεταγμένες της

list(pol3.geoms[1].interiors[0].coords)
[(30.0, 20.0), (20.0, 15.0), (20.0, 25.0), (30.0, 20.0)]

Ιδιότητες αντικειμένων

Υπολογισμός ορίων (bounds)

geo_collection
../_images/Lesson_9_124_0.svg
geo_collection.bounds
(-2.0, -1.0, 45.0, 45.0)
shapely.geometry.box(*geo_collection.bounds)
../_images/Lesson_9_126_0.svg
line1
../_images/Lesson_9_127_0.svg
line1.bounds
(-1.0, 0.0, 2.0, 1.0)
pnt1
../_images/Lesson_9_129_0.svg
pnt1.bounds
(2.0, 0.5, 2.0, 0.5)
list(pnt1.coords)
[(2.0, 0.5)]

Υπολογισμός μήκους γραμμής

line1
../_images/Lesson_9_133_0.svg
line1.length
5.354101966249685
line2.geoms[0]
../_images/Lesson_9_135_0.svg
line2.geoms[0].length
5.5901699437494745

Υπολογισμός εμβαδού

pol1.area
7.5
pol2.area
6.25

Νέες γεωμετρίες

Κεντροειδές πολυγώνου (centroid)

pol2
../_images/Lesson_9_142_0.svg
pol2.centroid
../_images/Lesson_9_143_0.svg
shapely.geometry.GeometryCollection([pol2, pol2.centroid])
../_images/Lesson_9_144_0.svg

Περιμετρική ζώνη (buffer)

pnt1.buffer(5)
../_images/Lesson_9_146_0.svg
shapely.geometry.GeometryCollection([pnt1,pnt1.buffer(5)])
../_images/Lesson_9_147_0.svg
pol1.buffer(5)
../_images/Lesson_9_148_0.svg
shapely.geometry.GeometryCollection([pol1,pol1.buffer(5)])
../_images/Lesson_9_149_0.svg

Convex hull

pol3
../_images/Lesson_9_151_0.svg
pol3.convex_hull
../_images/Lesson_9_152_0.svg

Σχέσεις μεταξύ αντικειμένων

shapely.geometry.GeometryCollection([pol1,pol3])
../_images/Lesson_9_154_0.svg
pol3
../_images/Lesson_9_155_0.svg
pol1.intersects(pol3)
False
pol1.intersects(pol2)
True
shapely.geometry.GeometryCollection([pol1,pol2])
../_images/Lesson_9_158_0.svg

Γεωμετρικές πράξεις

x = shapely.geometry.Point((0, 0)).buffer(1)
y = shapely.geometry.Point((1, 0)).buffer(1)
shapely.geometry.GeometryCollection([x, y])
../_images/Lesson_9_160_0.svg
x.intersection(y)
../_images/Lesson_9_161_0.svg
x.difference(y)
../_images/Lesson_9_162_0.svg
x.union(y)
../_images/Lesson_9_163_0.svg
pol1.union(pol2)
../_images/Lesson_9_164_0.svg

Υπολογισμός απόστασης ανάμεσα σε δύο αντικείμενα

shapely.geometry.GeometryCollection([pol1, pol3])
../_images/Lesson_9_166_0.svg
pol1.distance(pol3)
10.307764064044152
pol3.distance(pol1) 
10.307764064044152
import pyproj

from shapely.geometry import Point
from shapely.ops import transform

wgs84_pt = Point(39.35858398397631, 22.932070043920394)

wgs84 = pyproj.CRS('EPSG:4326')
greek_grid = pyproj.CRS('EPSG:2100')

project = pyproj.Transformer.from_crs(wgs84, greek_grid, always_xy=True).transform
projected_point = transform(project, wgs84_pt)
print(projected_point)
POINT (2087909.8290560723 2620022.621513317)
list(wgs84_pt.coords)[0]
(39.35858398397631, 22.932070043920394)
import folium

coords = list(wgs84_pt.coords)[0]

#Create the map
my_map = folium.Map(location = coords, zoom_start = 13)

# Add marker
folium.Marker(coords, popup = 'Volos').add_to(my_map)

#Display the map
my_map
Make this Notebook Trusted to load map: File -> Trust Notebook

Η βιβλιοθήκη Geopandas

Μέχρι τώρα είδαμε πως μπορούμε να διαχειριζόμαστε γεωμετρικά δεδομένα με την βιβλιοθήκη shapely.

Όμως τα γεωγραφικά δεδομένα δεν περιλαμβάνουν μόνο της γεωγραφική πληροφορία και την γεωμετρική τους δομή αλλά συνοδεύονται από μια σειρά περιγραφικών δεδομένων.

Η βιβλιοθήκη geopandas είναι μια βιβλιοθήκη της Python η οποία υποστηρίζει την ανάγνωση, επεξεργασία, ανάλυση και εγγραφή γεωγραφικών και παράλληλα περιγραφικών δεδομένων. Αποτελεί επέκταση της βιβλιοθήκης pandas και «κληρονομεί» τα χαρακτηριστικά και τις δυνατότητές της.

Στην υφιστάμενη δομή της pandas, η geopandas υποστηρίζει την γεωμετρία με την προσθήκη μιας νέας στήλης και το γεωγραφικό σύστημα αναφοράς.

alt text

alt text

Εισάγουμε την βιβλιοθήκη με τον παρακάτω τρόπο

import geopandas as gpd

Για να διαβάσουμε ένα αρχείο shapefile χρησιμοποιούμε την συνάρτηση gpd.read_file. Το αντικείμενο που προκύπτει είναι τύπου geopandas.geodataframe.GeoDataFrame

# Import shapefile using geopandas
dhmoi = gpd.read_file("../docs/dhmoi.gpkg")

type(dhmoi)
geopandas.geodataframe.GeoDataFrame

Η στήλη της γεωμετρίας ονομάζεται προκαθορισμένα geometry και είναι αντικείμενο geopandas.geoseries.GeoSeries που περιλαμβάνει αντικείμενα shapely.

type(dhmoi["geometry"])
geopandas.geoseries.GeoSeries

Με την μέθοδο geom_type μπορούμε να δούμε τον γεωμετρικό τύπο κάθε εγγραφής.

dhmoi.geom_type
0      MultiPolygon
1      MultiPolygon
2      MultiPolygon
3      MultiPolygon
4      MultiPolygon
           ...     
320    MultiPolygon
321    MultiPolygon
322    MultiPolygon
323    MultiPolygon
324    MultiPolygon
Length: 325, dtype: object

Ας δούμε τις αρχικές εγγραφές του αρχείου

dhmoi.head()
OBJECTID X Y Name CodeELSTAT PopM01 PopF01 PopTot01 UnemrM01 UnemrF01 UnemrT01 PrSect01 Foreig01 Income01 Perif geometry
0 1 616259.007833 4.551127e+06 KOMOTINIS 0101 29967.0 31534.0 61501.0 4.726948 6.051873 5.221281 21.817852 1.996780 10969.315256 1 MULTIPOLYGON (((631058.735 4574912.809, 631431...
1 2 644783.135624 4.561364e+06 ARRIANON 0102 8952.0 9307.0 18259.0 1.955813 2.390181 2.154240 92.912436 0.038337 6398.641486 1 MULTIPOLYGON (((661976.597 4554253.076, 662118...
2 3 602100.950461 4.555689e+06 IASMOU 0103 7290.0 7561.0 14851.0 5.653884 6.671554 6.062390 70.206767 1.339977 6608.897469 1 MULTIPOLYGON (((603623.489 4566376.502, 603657...
3 4 633712.468442 4.539548e+06 MARONEIAS - SAPON 0104 8284.0 8342.0 16626.0 5.840957 10.270499 7.421934 57.748085 0.890172 6700.375345 1 MULTIPOLYGON (((650998.525 4557469.860, 651446...
4 5 517831.195188 4.572544e+06 DRAMAS 0201 28041.0 29326.0 57367.0 8.831689 14.176926 10.851115 7.298838 2.961633 11299.106498 2 MULTIPOLYGON (((521986.700 4601855.354, 522271...
dhmoi.plot()
<AxesSubplot:>
../_images/Lesson_9_186_1.png

Χαρτογραφική απόδοση σε διαδραστικό χάρτη

dhmoi.explore(legend=False)
Make this Notebook Trusted to load map: File -> Trust Notebook

Εκτύπωση λεπτομερειών για το γεωγραφικό σύστημα αναφοράς το οποίο όπως φαίνεται παρακάτω είναι το ΕΓΣΑ “87 (EPSG:2100)

dhmoi.crs
<Derived Projected CRS: PROJCS["unknown",GEOGCS["unknown",DATUM["Greek_Geo ...>
Name: unknown
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Transverse Mercator
Datum: Greek Geodetic Reference System 1987
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

Μπορούμε να μετασχηματίσουμε τα δεδομένα σε ένα διαφορετικό προβολικό σύστημα με κλήση της μεθόδου to_crs

dhmoi_wgs84 = dhmoi.to_crs(4326)
dhmoi_wgs84.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
dhmoi.geom_type
0      MultiPolygon
1      MultiPolygon
2      MultiPolygon
3      MultiPolygon
4      MultiPolygon
           ...     
320    MultiPolygon
321    MultiPolygon
322    MultiPolygon
323    MultiPolygon
324    MultiPolygon
Length: 325, dtype: object

Εκτύπωση των γεωγραφικών ορίων του αρχείου oikismoi που είναι στο σύστημα αναφορά ΕΓΣΑ “87.

dhmoi.total_bounds
array([ 104040.7266, 3850938.    , 1007943.    , 4624010.2508])

Εκτύπωση των γεωγραφικών ορίων του αρχείου oikismoi που είναι στο σύστημα αναφορά WGS’84.

dhmoi_wgs84.total_bounds
array([19.37377196, 34.80060523, 29.64123782, 41.74911332])

Η μέθοδος shape μας επιστρέφει τον πλήθος των γραμμών (325) και των στηλών (15).

dhmoi.shape
(325, 16)

Μπορούμε να εγγράψουμε ένα Geopandas Dataframe

dhmoi.to_file("dhmoi.geojson", driver="GeoJSON")
/home/leonidas/anaconda3/envs/book/lib/python3.10/site-packages/geopandas/io/file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  pd.Int64Index,
dhmoi.sort_values(by='PopTot01', ascending=False)
OBJECTID X Y Name CodeELSTAT PopM01 PopF01 PopTot01 UnemrM01 UnemrF01 UnemrT01 PrSect01 Foreig01 Income01 Perif geometry
192 193 476676.870086 4.204338e+06 ATHINAION 4501 374900.0 414266.0 789166.0 4.890244 6.640842 5.624998 0.515649 17.426785 14497.393539 45 MULTIPOLYGON (((478164.675 4208572.203, 478015...
22 23 411331.939920 4.497210e+06 THESSALONIKIS 0701 185771.0 211385.0 397156.0 5.402357 7.615427 6.336333 0.834415 7.200951 13388.724605 7 MULTIPOLYGON (((408112.992 4500238.217, 408264...
147 148 306893.782874 4.231792e+06 PATREON 3701 104429.0 106065.0 210494.0 7.613916 9.985632 8.493895 2.567481 5.301864 12699.732926 37 MULTIPOLYGON (((312258.406 4245213.000, 312382...
245 246 468574.034133 4.199691e+06 PEIRAIOS 5101 87362.0 94571.0 181933.0 5.742561 7.683770 6.480184 0.491722 8.515141 13526.486840 51 MULTIPOLYGON (((469779.406 4198249.000, 469760...
301 302 599279.522964 3.900401e+06 IRAKLEIOU 7101 81067.0 82048.0 163115.0 5.047587 8.832661 6.557167 7.477305 4.048727 12479.911946 71 MULTIPOLYGON (((602945.188 3911414.000, 602942...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
299 300 824000.991323 4.016217e+06 CHALKIS 6905 165.0 130.0 295.0 1.265823 4.166667 1.941748 31.683168 17.966102 8888.138597 69 MULTIPOLYGON (((830115.688 4011995.000, 830198...
270 271 659987.832319 4.025277e+06 ANAFIS 6002 135.0 137.0 272.0 1.369863 2.941176 1.869159 43.809524 1.102941 9062.404166 60 MULTIPOLYGON (((669453.375 4014893.000, 669444...
272 273 599675.697252 4.059438e+06 SIKINOU 6004 118.0 120.0 238.0 0.000000 4.545455 1.333333 27.027027 1.680672 7566.666874 60 MULTIPOLYGON (((593612.000 4054459.000, 593596...
275 276 762387.857032 4.149989e+06 AGATHONISIOU 6102 92.0 60.0 152.0 12.765957 10.000000 12.280702 56.000000 9.210526 5421.733909 61 MULTIPOLYGON (((763781.625 4146430.000, 763859...
320 321 507562.964610 3.855768e+06 GAVDOU 7403 49.0 32.0 81.0 0.000000 0.000000 0.350000 51.351351 9.876543 8498.080328 74 MULTIPOLYGON (((505925.688 3858625.000, 506062...

325 rows × 16 columns

Ορισμός του index στην στήλη «CodeELSTAT»

dhmoi = dhmoi.set_index("CodeELSTAT")
dhmoi
OBJECTID X Y Name PopM01 PopF01 PopTot01 UnemrM01 UnemrF01 UnemrT01 PrSect01 Foreig01 Income01 Perif geometry
CodeELSTAT
0101 1 616259.007833 4.551127e+06 KOMOTINIS 29967.0 31534.0 61501.0 4.726948 6.051873 5.221281 21.817852 1.996780 10969.315256 1 MULTIPOLYGON (((631058.735 4574912.809, 631431...
0102 2 644783.135624 4.561364e+06 ARRIANON 8952.0 9307.0 18259.0 1.955813 2.390181 2.154240 92.912436 0.038337 6398.641486 1 MULTIPOLYGON (((661976.597 4554253.076, 662118...
0103 3 602100.950461 4.555689e+06 IASMOU 7290.0 7561.0 14851.0 5.653884 6.671554 6.062390 70.206767 1.339977 6608.897469 1 MULTIPOLYGON (((603623.489 4566376.502, 603657...
0104 4 633712.468442 4.539548e+06 MARONEIAS - SAPON 8284.0 8342.0 16626.0 5.840957 10.270499 7.421934 57.748085 0.890172 6700.375345 1 MULTIPOLYGON (((650998.525 4557469.860, 651446...
0201 5 517831.195188 4.572544e+06 DRAMAS 28041.0 29326.0 57367.0 8.831689 14.176926 10.851115 7.298838 2.961633 11299.106498 2 MULTIPOLYGON (((521986.700 4601855.354, 522271...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
7403 321 507562.964610 3.855768e+06 GAVDOU 49.0 32.0 81.0 0.000000 0.000000 0.350000 51.351351 9.876543 8498.080328 74 MULTIPOLYGON (((505925.688 3858625.000, 506062...
7404 322 475420.582085 3.905918e+06 KANTANOU - SELINOU 3271.0 3031.0 6302.0 2.312500 2.234637 2.284569 61.525841 14.680209 7989.958401 74 MULTIPOLYGON (((469969.812 3897335.000, 469955...
7405 323 465887.990746 3.920087e+06 KISSAMOU 5925.0 5545.0 11470.0 4.847458 5.894591 5.191257 45.701249 7.916303 8750.030129 74 MULTIPOLYGON (((457703.500 3902963.000, 457852...
7406 324 484468.969583 3.924051e+06 PLATANIA 9217.0 8647.0 17864.0 4.437401 8.586262 5.938448 51.305684 8.889386 8120.027122 74 MULTIPOLYGON (((476162.812 3949684.000, 476253...
7407 325 508083.904768 3.901847e+06 SFAKION 1288.0 1131.0 2419.0 1.920439 2.839117 2.198853 63.929619 4.960728 6602.391386 74 MULTIPOLYGON (((502587.434 3912125.553, 502774...

325 rows × 15 columns

Υπολογισμός έκτασης (σε \(km^2\)) σε μία νέα στήλη με το όνομα area

dhmoi["area"] = dhmoi.area*10e-6
dhmoi
OBJECTID X Y Name PopM01 PopF01 PopTot01 UnemrM01 UnemrF01 UnemrT01 PrSect01 Foreig01 Income01 Perif geometry area
CodeELSTAT
0101 1 616259.007833 4.551127e+06 KOMOTINIS 29967.0 31534.0 61501.0 4.726948 6.051873 5.221281 21.817852 1.996780 10969.315256 1 MULTIPOLYGON (((631058.735 4574912.809, 631431... 6470.182562
0102 2 644783.135624 4.561364e+06 ARRIANON 8952.0 9307.0 18259.0 1.955813 2.390181 2.154240 92.912436 0.038337 6398.641486 1 MULTIPOLYGON (((661976.597 4554253.076, 662118... 7724.393191
0103 3 602100.950461 4.555689e+06 IASMOU 7290.0 7561.0 14851.0 5.653884 6.671554 6.062390 70.206767 1.339977 6608.897469 1 MULTIPOLYGON (((603623.489 4566376.502, 603657... 4870.493148
0104 4 633712.468442 4.539548e+06 MARONEIAS - SAPON 8284.0 8342.0 16626.0 5.840957 10.270499 7.421934 57.748085 0.890172 6700.375345 1 MULTIPOLYGON (((650998.525 4557469.860, 651446... 6441.606153
0201 5 517831.195188 4.572544e+06 DRAMAS 28041.0 29326.0 57367.0 8.831689 14.176926 10.851115 7.298838 2.961633 11299.106498 2 MULTIPOLYGON (((521986.700 4601855.354, 522271... 8401.321201
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
7403 321 507562.964610 3.855768e+06 GAVDOU 49.0 32.0 81.0 0.000000 0.000000 0.350000 51.351351 9.876543 8498.080328 74 MULTIPOLYGON (((505925.688 3858625.000, 506062... 343.773091
7404 322 475420.582085 3.905918e+06 KANTANOU - SELINOU 3271.0 3031.0 6302.0 2.312500 2.234637 2.284569 61.525841 14.680209 7989.958401 74 MULTIPOLYGON (((469969.812 3897335.000, 469955... 3746.368061
7405 323 465887.990746 3.920087e+06 KISSAMOU 5925.0 5545.0 11470.0 4.847458 5.894591 5.191257 45.701249 7.916303 8750.030129 74 MULTIPOLYGON (((457703.500 3902963.000, 457852... 3421.398731
7406 324 484468.969583 3.924051e+06 PLATANIA 9217.0 8647.0 17864.0 4.437401 8.586262 5.938448 51.305684 8.889386 8120.027122 74 MULTIPOLYGON (((476162.812 3949684.000, 476253... 4925.063676
7407 325 508083.904768 3.901847e+06 SFAKION 1288.0 1131.0 2419.0 1.920439 2.839117 2.198853 63.929619 4.960728 6602.391386 74 MULTIPOLYGON (((502587.434 3912125.553, 502774... 4667.779146

325 rows × 16 columns

Υπολογισμός του centroid κάθε δήμου σε μια νέα στήλη (centroid)

dhmoi['centroid'] = dhmoi.centroid
dhmoi
OBJECTID X Y Name PopM01 PopF01 PopTot01 UnemrM01 UnemrF01 UnemrT01 PrSect01 Foreig01 Income01 Perif geometry area centroid
CodeELSTAT
0101 1 616259.007833 4.551127e+06 KOMOTINIS 29967.0 31534.0 61501.0 4.726948 6.051873 5.221281 21.817852 1.996780 10969.315256 1 MULTIPOLYGON (((631058.735 4574912.809, 631431... 6470.182562 POINT (616259.008 4551127.241)
0102 2 644783.135624 4.561364e+06 ARRIANON 8952.0 9307.0 18259.0 1.955813 2.390181 2.154240 92.912436 0.038337 6398.641486 1 MULTIPOLYGON (((661976.597 4554253.076, 662118... 7724.393191 POINT (644783.136 4561364.036)
0103 3 602100.950461 4.555689e+06 IASMOU 7290.0 7561.0 14851.0 5.653884 6.671554 6.062390 70.206767 1.339977 6608.897469 1 MULTIPOLYGON (((603623.489 4566376.502, 603657... 4870.493148 POINT (602100.950 4555689.138)
0104 4 633712.468442 4.539548e+06 MARONEIAS - SAPON 8284.0 8342.0 16626.0 5.840957 10.270499 7.421934 57.748085 0.890172 6700.375345 1 MULTIPOLYGON (((650998.525 4557469.860, 651446... 6441.606153 POINT (633712.468 4539548.490)
0201 5 517831.195188 4.572544e+06 DRAMAS 28041.0 29326.0 57367.0 8.831689 14.176926 10.851115 7.298838 2.961633 11299.106498 2 MULTIPOLYGON (((521986.700 4601855.354, 522271... 8401.321201 POINT (517831.195 4572544.478)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
7403 321 507562.964610 3.855768e+06 GAVDOU 49.0 32.0 81.0 0.000000 0.000000 0.350000 51.351351 9.876543 8498.080328 74 MULTIPOLYGON (((505925.688 3858625.000, 506062... 343.773091 POINT (507562.965 3855767.657)
7404 322 475420.582085 3.905918e+06 KANTANOU - SELINOU 3271.0 3031.0 6302.0 2.312500 2.234637 2.284569 61.525841 14.680209 7989.958401 74 MULTIPOLYGON (((469969.812 3897335.000, 469955... 3746.368061 POINT (475420.582 3905917.949)
7405 323 465887.990746 3.920087e+06 KISSAMOU 5925.0 5545.0 11470.0 4.847458 5.894591 5.191257 45.701249 7.916303 8750.030129 74 MULTIPOLYGON (((457703.500 3902963.000, 457852... 3421.398731 POINT (465887.991 3920086.947)
7406 324 484468.969583 3.924051e+06 PLATANIA 9217.0 8647.0 17864.0 4.437401 8.586262 5.938448 51.305684 8.889386 8120.027122 74 MULTIPOLYGON (((476162.812 3949684.000, 476253... 4925.063676 POINT (484468.970 3924050.753)
7407 325 508083.904768 3.901847e+06 SFAKION 1288.0 1131.0 2419.0 1.920439 2.839117 2.198853 63.929619 4.960728 6602.391386 74 MULTIPOLYGON (((502587.434 3912125.553, 502774... 4667.779146 POINT (508083.905 3901847.180)

325 rows × 17 columns

Χαρτογραφική απόδοση του δείκτη ανεργίας (στήλη UnemrT01) ανά δήμο

dhmoi.plot("UnemrT01", legend=True)
<AxesSubplot:>
../_images/Lesson_9_213_1.png
 dhmoi.plot("Income01", scheme='quantiles', cmap='YlOrRd',  legend=True, figsize=(12, 10))
<AxesSubplot:>
../_images/Lesson_9_214_1.png

Μπορούμε να οπτικοποιήσουμε διαφορετική στήλη τύπου geopandas.geoseries.GeoSeries αφού πρώτα την ορίσουμε με την μέθοδο set_geometry.

dhmoi = dhmoi.set_geometry("centroid")
dhmoi.plot("UnemrT01", legend=True)
<AxesSubplot:>
../_images/Lesson_9_216_1.png

Ορισμός περιμετρικής ζώνης διαμέτρου 20χλμ γύρω από κάθε centroid.

dhmoi = dhmoi.set_geometry("centroid") # ορισμός default γεωμετρίας η στήλη centroid
dhmoi["buffered"] = dhmoi.buffer(20000) #εκτέλεση περιμετρικών ζωνών

dhmoi = dhmoi.set_geometry("buffered") # ορισμός default γεωμετρίας η στήλη buffered
dhmoi.plot(legend=True) # οπτικοποίηση
<AxesSubplot:>
../_images/Lesson_9_218_1.png

Ξανα ορίζουμε σαν προκαθορισμένη στήλη γεωμετρία την στήλη geometry.

dhmoi = dhmoi.set_geometry("geometry")

Μπορούμε να φιλτράρουμε γραμμές και στήλες όπως και στην pandas χρησιμοποιώντας το index και το όνομα στήλης:

lesvos = dhmoi.loc["5301", "geometry"]
lesvos
../_images/Lesson_9_223_0.svg
type(lesvos)
shapely.geometry.multipolygon.MultiPolygon

Η παραπάνω διαδικασία επιλογής επέστρεψε ένα shapely MultiPolygon object όπως φαίνεται.

Με την μέθοδο dissolve μπορούμε να συγχωνεύσουμε γεωμετρικά αντικείμενα βάση κάποιας στήλης. Αντικείμενα με ίδια τιμή θα συγχωνευθούν σε ενιαίο γεωμετρικό αντικείμενο.

dhmoi.head()
OBJECTID X Y Name PopM01 PopF01 PopTot01 UnemrM01 UnemrF01 UnemrT01 PrSect01 Foreig01 Income01 Perif geometry area centroid buffered
CodeELSTAT
0101 1 616259.007833 4.551127e+06 KOMOTINIS 29967.0 31534.0 61501.0 4.726948 6.051873 5.221281 21.817852 1.996780 10969.315256 1 MULTIPOLYGON (((631058.735 4574912.809, 631431... 6470.182562 POINT (616259.008 4551127.241) POLYGON ((636259.008 4551127.241, 636162.702 4...
0102 2 644783.135624 4.561364e+06 ARRIANON 8952.0 9307.0 18259.0 1.955813 2.390181 2.154240 92.912436 0.038337 6398.641486 1 MULTIPOLYGON (((661976.597 4554253.076, 662118... 7724.393191 POINT (644783.136 4561364.036) POLYGON ((664783.136 4561364.036, 664686.830 4...
0103 3 602100.950461 4.555689e+06 IASMOU 7290.0 7561.0 14851.0 5.653884 6.671554 6.062390 70.206767 1.339977 6608.897469 1 MULTIPOLYGON (((603623.489 4566376.502, 603657... 4870.493148 POINT (602100.950 4555689.138) POLYGON ((622100.950 4555689.138, 622004.645 4...
0104 4 633712.468442 4.539548e+06 MARONEIAS - SAPON 8284.0 8342.0 16626.0 5.840957 10.270499 7.421934 57.748085 0.890172 6700.375345 1 MULTIPOLYGON (((650998.525 4557469.860, 651446... 6441.606153 POINT (633712.468 4539548.490) POLYGON ((653712.468 4539548.490, 653616.163 4...
0201 5 517831.195188 4.572544e+06 DRAMAS 28041.0 29326.0 57367.0 8.831689 14.176926 10.851115 7.298838 2.961633 11299.106498 2 MULTIPOLYGON (((521986.700 4601855.354, 522271... 8401.321201 POINT (517831.195 4572544.478) POLYGON ((537831.195 4572544.478, 537734.890 4...
nomoi = dhmoi.dissolve(by='Perif')
dhmoi.plot()
<AxesSubplot:>
../_images/Lesson_9_229_1.png
nomoi.plot()
<AxesSubplot:>
../_images/Lesson_9_230_1.png

Επιπλεόν κατά την συγχώνευση μπορούμε να εφαρμόσουμε μια συνάρτηση στις στήλες π.χ. mean ή sum.

nomoi = dhmoi.dissolve(by='Perif', aggfunc='mean')
nomoi.plot(column = 'Income01', scheme='quantiles', cmap='YlOrRd');
../_images/Lesson_9_233_0.png

Βιβλιογραφία: